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Abstract 

In this paper we present a method to treat interface jump conditions for 
constant coefficients Poisson problems that allows the use of standard "black 
box" solvers, without compromising accuracy. The basic idea of the new ap- 
proach is similar to the Ghost Fluid Method (GFM). The GFM relies on cor- 
rections applied on nodes located across the interface for discretization sten- 
cils that straddle the interface. If the corrections are solution- independent, 
they can be moved to the right-hand-side (RHS) of the equations, producing 
a problem with the same linear system as if there were no jumps, only with a 
different RHS. However, achieving high accuracy is very hard (if not impos- 
sible) with the "standard" approaches used to compute the GFM correction 
terms. 

In this paper we generalize the GFM correction terms to a correction func- 
tion, defined on a band around the interface. This function is then shown to 
be characterized as the solution to a PDE, with appropriate boundary condi- 
tions. This PDE can, in principle, be solved to any desired order of accuracy. 
As an example, we apply this new method to devise a 4**^ order accurate 
scheme for the constant coefficients Poisson equation with discontinuities in 
2D. This scheme is based on (i) the standard 9-point stencil discretization 
of the Poisson equation, (ii) a representation of the correction function in 
terms of bicubics, and (iii) a solution of the correction function PDE by a 
least squares minimization. Several applications of the method are presented 
to illustrate its robustness dealing with a variety of interface geometries, its 
capability to capture sharp discontinuities, and its high convergence rate. 
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1. Introduction. 

1.1. Motivation and background information. 

In this paper we present a new and efficient method to solve the constant 
coefficients Poisson equation in the presence of discontinuities across an in- 
terface, with a high order of accuracy. Solutions of the Poisson equation with 
discontinuities are of fundamental importance in the description of fluid flows 
separated by interfaces {e.g. the contact surfaces for immiscible multiphase 
fluids, or fluids separated by a membrane) and other multiphase diffusion 
phenomena. Over the last three decades, several methods have been devel- 
oped to solve problems of this type numerically p14T7]. However, obtaining 
a high order of accuracy still poses great challenges in terms of complexity 
and computational efficiency. 

When the solution is known to be smooth, it is easy to obtain highly 
accurate finite-difference discretizations of the Poisson equation on a regu- 
lar grid. Furthermore, these discretizations commonly yield symmetric and 
banded linear systems, which can be inverted efficiently [18]. On the other 
hand, when singularities occur {e.g. discontinuities) across internal interfaces, 
some of the regular discretization stencils will straddle the interface, which 
renders the whole procedure invalid. 

Several strategies have been proposed to tackle this issue. Peskin [T] intro- 
duced the Immersed Boundary Method (IBM) [H , in which the disconti- 
nuities are re-interpreted as additional (singular) source terms concentrated 
on the interface. These singular terms are then "regularized" and appro- 
priately spread out over the regular grid — in a "thin" band enclosing the 
interface. The result is a ffist order scheme that smears discontinuities. In 
order to avoid this smearing of the interface information, LeVeque and Li [3] 
developed the Immersed Interface Method (IIM) [3l SI [HI US], which is a 
methodology to modify the discretization stencils, taking into consideration 
the discontinuities at their actual locations. The IIM guarantees second order 
accuracy and sharp discontinuities, but at the cost of added discretization 
complexity and loss of symmetry. 
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The new method advanced in this paper builds on the ideas introduced 
by the Ghost Fluid Method (GFM) [SHSl [El HI HE]. The GFM is based 
on defining both actual and "ghost" fiuid variables at every node on a nar- 
row band enclosing the interface. The ghost variables work as extensions of 
the actual variables across the interface — the solution on each side of the 
interface is assumed to have a smooth extension into the other side. This ap- 
proach allows the use of standard discretizations everywhere in the domain. 
In most GFM versions, the ghost values are written as the actual values, plus 
corrections that are independent of the underlying solution to the Poisson 
problem. Hence, the corrections can be pre-computed, and moved into the 
source term for the equation. In this fashion the GFM yields the same linear 
system as the one produced by the problem without an interface, except for 
changes in the right-hand-side (sources) only, which can then be inverted just 
as efficiently. 

The key difficulty in the GFM is the calculation of the correction terms, 
since the overall accuracy of the scheme depends heavily on the quality of the 
assigned ghost values. In [BHOl [12] the authors develop first order accurate 
approaches to deal with discontinuities. In the present work, we show that 
for the constant coefficients Poisson equation we can generalize the GFM cor- 
rection term (at each ghost point) concept to that of a correction function 
defined on a narrow band enclosing the interface. Hence we call this new 
approach the Correction Function Method (GFM). This correction function 
is then shown to be characterized as the solution to a PDE, with appro- 
priate boundary conditions on the interface — see § |4j Thus, at least in 
principle, one can calculate the correction function to any order of accuracy, 
by designing algorithms to solve the PDE that defines it. In this paper we 
present examples of 2'^'^ and 4**^ order accurate schemes (to solve the con- 
stant coefficients Poisson equation, with discontinuities across interfaces, in 
2D) developed using this general framework. 

A key point (see § [s]) in the scheme developed here is the way we solve 
the PDE defining the correction function. This PDE is solved in a weak 
fashion using a least squares minimization procedure. This provides a flexible 
approach that allows the development of a robust scheme that can deal with 
the geometrical complications of the placement of the regular grid stencils 
relative to the interface. Furthermore, this approach is easy to generalize to 
3D, or to higher orders of accuracy. 
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1.2. Other related work. 

It is relevant to note other developments to solve the Poisson equation 
under similar circumstances — multiple phases separated by interfaces — 
but with different interface conditions. The Poisson problem with Dirichlet 
boundary conditions, on an irregular boundary embedded in a regular grid, 
has been solved to second order of accuracy using fast Poisson solvers [T9] . 
finite volume [3], and finite differences [20H23] approaches. In particular, 
Gibou et al. [21] and Jomaa and Macaskill |22] have shown that it is possi- 
ble to obtain symmetric discretizations of the embedded Dirichlet problem, 
up to second order of accuracy. Gibou and Fedkiw [23] have developed a 
fourth order accurate discretization of the problem, at the cost of giving up 
symmetry. More recently, the same problem has also been solved, to second 
order of accuracy, in non-graded adaptive Cartesian grids by Chen et al. . 
Furthermore, the embedded Dirichlet problem is closely related to the Stefan 
problem modeling dendritic growth, as described in [251 126] . 

The finite-element community has also made significant progress in incor- 
porating the IIM and similar techniques to solve the Poisson equation using 
embedded grids. In particular the works by Gong et al. [T3], Dolbow and 
Harari [T6|, and Bedrossian et al. [T7] describe second order accurate finite- 
element discretizations that result in symmetric linear systems. Moreover, 
in these works the interface (or boundary) conditions are imposed in a weak 
fashion, which bears some conceptual similarities with the CFM presented 
here, although the execution is rather different. 

1.3. Interface representation. 

Another issue of primary importance to multiphase problems is the repre- 
sentation of the interface (and its tracking in unsteady cases). Some authors 
(see [H [ini EH]) choose to represent the interface explicitly, by tracking in- 
terface particles. The location of the neighboring particles is then used to 
produce local interpolations {e.g. splines), which are then applied to compute 
geometric information — such as curvature and normal directions. Although 
this approach can be quite accurate, it requires special treatment when the 
interface undergoes either large deformations or topological changes — such 
as mergers or splits. Even though we are not concerned with these issues 
in this paper, we elected to adopt an implicit representation, to avoid com- 
plications in future applications. In an implicit representation, the interface 
is given as the zero level of a function that is defined everywhere in the 
regular grid — the level set function [27]. In particular, we adopted the 
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Gradient-Augmented Level Set (GA-LS) method [2B]- With this extension 
of the level set method, we can obtain highly accurate representations of the 
interface, and other geometric information, with the additional advantage 
that this method uses only local grid information. We discuss the question 



of interface representation in more detail in § 5.6 



1.4- Organization of the paper. 

The remainder of the paper is organized as follows. In § [2] we introduce 
the Poisson problem that we seek to solve. In § [3] the basic idea behind 
the solution method, and its relationship to the GFM, are explored. Next, 
in § |4| we introduce the concept of the correction function and show how 
it is defined by a PDE problem. In § [5] we apply this new framework to 
build a 4*"^ order accurate scheme in 2D. Since the emphasis of this paper 
is on high-order schemes, we describe the 2"*^ order accurate scheme in ap- 
pendix [Cj Next, in § [6] we demonstrate the robustness and accuracy of the 
2D scheme by applying it to several example problems. The conclusions are 
in § [7j In appendix |A] we review some background material, and notation, 
on bicubic interpolation. Finally, in appendix [B] we discuss some technical 
issues regarding the construction of the sets where the correction function is 
solved for. 



2. Definition of the problem. 

Our objective is to solve the constant coefficients Poisson's equation in 
a domain Q in which the solution is discontinuous across a co-dimension 
1 interface F, which divides the domain into the subdomains and Q~ , 
as illustrated in figure [H We use the notation and u to denote the 
solution in each of the subdomains. Let the discontinuities across F be given 
in terms of two functions defined on the interface: a = a{x) for the jump 
in the function values, and b = b{x) for the jump in the normal derivatives. 
Furthermore, assume Dirichlet boundary conditions on the "outer" boundary 
dVl (see figure [T]). Thus the problem to be solved is 

V^M (f) = / (f) for X eQ, (la) 

[u]t = a (x) for x eT, (lb) 

= b (x) for X eT, (Ic) 

u{x) = g (x) for x G dVL, (Id) 



5 



where 

[u]r = (x) — (x) for x G F, (2a) 

[u„]p = {%) — (x) for a; G r. (2b) 

Throughout this paper, x = (xi, 0:2, . . . ) ^ is the spatial vector (where 
z/ = 2, or = 3), and is the Laplacian operator defined by 

- E Ij- (3) 

i=l ' 

Furthermore, 

Un = n-Vu = h- {u^^, u^^, ...) (4) 

denotes the derivative of u in the direction of n, the unit vector normal to 
the interface T pointing towards (see figure 




Figure 1: Example of a domain f2, with an interface T. 

It is important to note that our method focuses on the discretization 
of the problem in the vicinity of the interface only. Thus, the method is 
compatible with any set of boundary conditions on dQ, not just Dirichlet. 
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3. Solution method — the basic idea. 

To achieve the goal of efficient high order discretizations of Poisson equa- 
tion in the presence of discontinuities, we build on the idea of the Ghost Fluid 
Method (GFM). In essence, we use a standard discretization of the Laplace 
operator (on a domain without an interface F) and modify the right-hand- 
side (RHS) to incorporate the jump conditions across F. Thus, the resulting 
linear system can be inverted as efficiently as in the case of a solution without 
discontinuities. 




Figure 2: Example in ID of a solution with a jump discontinuity. 

Let US ffist illustrate the key concept in the GFM with a simple example, 
involving a regular grid and the standard second order discretization of the 
ID analog of the problem we are interested in. Thus, consider the problem of 
discretizing the equation Uxx = f{x) in some interval Q = {x : xl < x < xr}, 
where u is discontinuous across some point xr — hence (respectively Q~) 
is the domain xl < x < xr (respectively xr < x < xr). Then, see figure |2} 
when trying to approximate u^x at a grid point Xj such that Xi < xy < a^^+i, 
we would like to write 

+ _ n+ 1 - 2uj + -»+ ^ 

where h the grid spacing. However, we do not have information 

on M^i, but rather on Thus, the idea is to estimate a correction for 
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to recover it^x? such that Eq. ([5]) can be apphed: 



, ^ <i-2n.+ +(n,^i + A+i) 

where -Dj+i = m^j^ — is the correction term. Now we note that, if -Dj+i 
can be written as a correction that is independent on the solution u, then it 
can be moved to the RHS of the equation, and absorbed into /. That is 

This allows the solution of the problem with prescribed discontinuities using 
the same discretization as the one employed to solve the simpler problem 
without an interface — which leads to a great efficiency gain. 

Remark 1. The error in estimating D is crucial in determining the accuracy 
of the final discretization. Liu, Fedkiw, and Kang [9J introduced a dimension- 
by-dimension linear extrapolation of the interface jump conditions, to get a 
first order approximation for D. Our new method is based on generalizing 
the idea of a correction term to that of a correction function, for which we 
can write an equation. One can then obtain high accuracy representations for 
D by solving this equation, without the complications into which dimension- 
by-dimension (with Taylor expansions) approaches run into. X 

Remark 2. An additional advantage of the correction function approach is 
that D can be calculated at any point near the interface F. Hence it can 
be used with any finite differences discretization of the Poisson equation, 
without regard to the particulars of its stencil (as would be the case with 
any approach based on Taylor expansions). 4> 



4. The correction function and the equation defining it. 

As mentioned earlier, the aim here is to generalize the correction term 
concept to that of a correction function, and then to find an equation (a 
PDE, with appropriate boundary conditions) that uniquely characterizes the 
correction function. Then, at least in principle, one can design algorithms to 
solve the PDE in order to obtain solutions to the correction function of any 
desired order of accuracy. 
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Let us begin by considering a small region Q-p enclosing the interface F, 
defined as the set of all the points within some distance TZ of F, where TZ is 
of the order of the grid size h. As we will see below, we would like to have 
TZ as small as possible. On the other hand, fir has to include all the points 
where the CFM requires corrections to be computed, which mean that 7^ 
cannot be smaller than \/2h. In addition, algorithmic considerations (to be 
seen later) force TZ to be slightly larger than this last value. 

Next, we assume we that can extrapolate both and u~ , so that they 
are valid everywhere within Qr, in such a way that they satisfy the Poisson 
equations 

V\+ (f) = /+ (x) for xe^r, (8a) 

(x) = f- (x) for X G fir, (8b) 

where and /~ are smooth enough (see remark |3] below) extensions of the 
source term / to fir- In particular, notice that the introduction of and 
/~ allows the possibility of the source term changing [i.e. a discontinuous 
source term) across F. The correction function is then defined by D{x) = 
u~^{x) — u~{x). 

Taking the difference between the equations in dsj), and using the jump 



conditions (lb Ic), yields 



V^D (x) = /+ (x) - f- (x) = Jd (x) for f G fir, (9a) 

D{x) =a (x) for f G F, (9b) 

Dn (x) = b (x) for X G F. (9c) 

This achieves the aim of having the correction function defined by a set of 
equations, with some provisos — see remark |4] below. Note that: 

1. If /+ (x) = /" (x), for X G fir, then (x) = 0, for x G fir- 

2. Equation ( [9c) ) imposes the true jump condition in the normal direction, 
whereas some versions of the GFM rely on a dimension-by-dimension 
approximation of this condition (see Ref. [H]). 

Remark 3. The smoothness requirement on and /~ is tied up to how 
accurate an approximation to the correction term D is needed. For example. 



For the particular discretization of the Laplace operator that we use in this paper. 
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if a 4**^ order algorithm is used to find D, this will (generally) necessitate D 
to be at least for the errors to actually be 4**^ order. Hence, in this case, 
fD = f+- f- must be C^. X 

Remark 4. Equation ^ is an elliptic Cauchy problem for D in ilr- In 
general, such problems are ill-posed. However, we are seeking for solutions 
within the context of a numerical approximation where 

(a) There is a frequency cut-off in both the data a = a{x) and b = b{x), and 
the description of the curve F. 

(b) We are interested in the solution only a small distance away from the 
interface F, where this distance vanishes simultaneously with the inverse 
of the cut-off frequency in point (a). 

What (a) and (b) mean is that the arbitrarily large growth rate for arbitrarily 
small perturbations, which is responsible for the ill-posedness of the Cauchy 
problem in Eq. (|9]), does not occur within the special context where we need 
to solve the problem. This large growth rate does not occur because, for the 
solutions of the Poisson equation, the growth rate for a perturbation of wave 
number < < oo along some straight line, is given by e'^'"^'^ — where d is 
the distance from the line. However, by construction, in the case of interest 
to us kd is bounded. Jit 

Remark 5. Let us be more precise, and define a number characterizing how 
well posed the discretized version of Eq. Q is, by 

a = largest growth rate possible, 

where growth is defined relative to the size of a perturbation to the solution 
on the interface. This number is determined by TZ (the "radius" of Qr) as 
the following calculation shows: First of all, there is no loss of generality in 
assuming that the interface is fiat, provided that the numerical grid is fine 
enough to resolve F. In this case, let us introduce an orthogonal coordinate 
system y on F, and let d be the signed distance to F (say, d > in Q~). 
Expanding the perturbations in Fourier modes along the interface, the typical 
mode has the form 

2TTik-y±2nkd 

— » — # 

where k is the Fourier wave vector, and k = \k\. The shortest wave-length 
that can be represented on a grid with mesh size < h <^ 1 corresponds to 
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k = fcmax = 1/ (2/i). Hence, we obtain the estimate 

for the maximum growth rate. X 

Remark 6. Clearly, a is intimately related to the condition number for the 
discretized problem — see § [5] In fact, at leading order, the two numbers 
should be (roughly) proportional to each other — with a proportionality con- 
stant that depends on the details of the discretization. For the discretization 
used in this paper (described further below), \/2h <TZ< 2\/2h, which leads 
to the rough estimate 85 < a < 7, 200. On the other hand, the observed con- 
dition numbers vary between 5,000 and 10,000. Hence, the actual condition 
numbers are only slightly higher than a for the ranges of grid sizes h that 
we used (we did not explore the asymptotic limit h 0). X 

Remark 7. Equation ^ depends on the known inputs for the problem only. 
Namely: /"*", /~, a, and b. Consequently D does not depend on the solution 
u. Hence, after solving for D, we can use a discretization for u that does not 
involve the interface: Whenever u is discontinuous, we evaluate D where the 
correction is needed, and transfer these values to the RHS. 4> 



Remark 8. When developing an algorithm for a linear Cauchy problem, 
such as the one in Eq. ([9]), the two key requirements are consistency and 
stability. In particular, when the solution depends on the "initial conditions" 
globally, stability (typically) imposes stringent constraints on the "time" step 
for any local (explicit) scheme. This would seem to suggest that, in order 
to solve Eq. ([9]), a "global" (involving the whole domain f2r) method will 
be needed. This, however, is not true: because we need to solve Eq. (|9| for 
one "time" step only — i.e. within an 0{h) distance from F, stability is not 
relevant. Hence, consistency is enough, and a fully local scheme is possible. 
In the algorithm described in § [s] we found that, for (local) quadrangular 
patches, the Cauchy problem leads to a well behaved algorithm when the 
length of the interface contained in each patch is of the same order as the 
diagonal length of the patch. This result is in line with the calculation in 
remark |5| we want to keep the "wavelength" (along F) of the perturbations 
introduced by the discretization as long as possible. In particular, this should 
then minimize the condition number for the local problems — see remark|6} X 
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5. A 4*"^ Order Accurate Scheme in 2D. 
5.1. Overview. 

In this section we use the general ideas presented earlier to develop a 
specific example of a 4*^ order accurate scheme in 2D. Before proceeding 
with an in-depth description of the scheme, we highlight a few key points: 

(a) We discretize Poisson's equation using a compact 9-point stencil. Com- 
pactness is important since it is directly related to the size of TZc, which 
has a direct impact on the problem's conditioning — see remarks |4] - [6] 

(b) We approximate D using bicubic interpolations (bicubics), each valid in a 
small neighborhood fip-' of the interface. This guarantees local 4**^ order 
accuracy with only 12 interpolation parameters — see [2H]- Each Qp-' 
corresponds to a point in the grid at which the standard discretization 
of Poisson's equation involves a stencil that straddles the interface P. 

(c) The domains f2p-' are rectangular regions, each enclosing a portion of P, 
and all the nodes where D is needed to complete the discretization of the 
Poisson equation at the (z, j)-th stencil. Each is a sub-domain of Qr- 

(d) Starting from (b) and (c), we design a local solver that provides an 
approximation to D inside each domain f2p"'. 

(e) The interface P is represented using the Gradient-Augmented Level Set 
approach — see |28] . This guarantees a local 4*^ order representation of 
the interface, as required to keep the overall accuracy of the scheme. 

(f) In each Q^^ , we solve the PDE in Q in a least squares sense. Namely: 
First we define an appropriate positive quadratic integral quantity Jp — 



Eq. (17) — for which the solution is a minimum (actually, zero). Next 
we substitute the bicubic approximation for the solution into Jp, and 
discretize the integrals using Gaussian quadrature. Finally, we find the 
bicubic parameters by minimizing the discretized Jp. 

Remark 9. Solving the PDE in a least squares sense is crucial, since an 
algorithm is needed that can deal with the myriad ways in which the inter- 
face P can be placed relative to the fixed rectangular grid used to discretize 
Poisson's equation. This approach provides a scheme that (i) is robust with 
respect to the details of the interface geometry, (ii) has a formulation that 
is (essentially) dimension independent — there are no fundamental changes 
from 2D to 3D, and (iii) has a clear theoretical underpinning that allows 
extensions to higher orders, or to other discretizations of the Poisson equa- 
tion. 4> 



12 



5.2. Standard Stencil. 

We use the standard 4**^ order accurate 9-point discretization of Poisson's 
equatior]|^ 



L^'^ij + ^ (^X + hi) d^xdyyUi 



+ ^ {hi {f..kj + hi Uyy)i,) , (10) 



where is the second order 5-point discretization of the Laplace operator: 



and 



9xx^i,j 



'^yy'^hj 



h^x ' 



(11) 

(12) 
(13) 



The terms {fxx)i,j and {fyy)ij may be given analytically (if known), or com- 
puted using appropriate second order discretizations. 

In the absence of discontinuities, Eq. (10) provides a compact 4*"^ order 



accurate representation of Poisson's equation. In the vicinity of the disconti- 
nuities at the interface F, we define an appropriate domain •> ^^"^ compute 
the correction terms necessary to Eq. (10) — as described in detail next. 



To understand how the correction terms affect the discretization, let us 
consider the situation depicted in figure [3] In this case, the node (z,j) lies in 

while the nodes {i + 1, j), {i + 1, j + 1), and (z, j + 1) are in VL~ . Hence, 
to be able to use Eq. (10), we need to compute -Di+ij, -Dj+ij+i, and -Dij+i- 



After having solved for D where necessary (see § 5.3 and § 5.4), we modify 



Eq. ( 10 ) and write 



L^'uij + ^{hl + hi) dxxd, 



yy'^ij 



fij + {fxx)i,j + hy {fyy)i,j) + Cjj, 

(14) 

which differs from Eq. (10) by the terms Cij on the RHS only. Here the 



Cij are the CFM correction terms needed to complete the stencil across the 



^Notice that here we aUow for the possibihty of different grid spacings in each direction. 
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Figure 3: The 9-point compact stencil next to the interface T. 

discontinuity at T. In the particular case illustrated in figure [3| we have 
_\l {hl + hl) 1 

1 (hi + hi) 
12 ih^hyf ""'^''^ 

Similar formulas apply for the other possible arrangements of the Poisson's 
equation stencil relative to the interface F. 

5.3. Definition o/fip-'. 

There is some freedom on how to define fip"'. The basic requirements are 

(i) fl^Y should be a rectangle. 

(ii) the edges of fl^^'' should be parallel to the grid lines. 

(iii) Qp-' should be small, since the problem's condition number increases 
exponentially with the distance from F — see remarks [5] and [6j 

(iv) fip"* should contain all the nodes where D is needed. For the example, 
in figure |3] we need to know Dj+ij+i, A+ij, and Aj+i- Hence, in this 
case, f2p"' should include the nodes (i + l,j + l), (i + l,j), and (i,j + l). 



1 {hi + hi) 
6 



{h^h 



y) 



hi 



D 
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(v) fip-' should contain a segment of F, with a length that is as large as 
possible — i.e. comparable to the length of the diagonal of Q^'' ■ This 
follows from the calculation in remark |5| which indicates that the wave- 
length of the perturbations (along F) introduced by the discretization 
should be as long as possible. This should then minimize the condition 
number for the local problem — see remark [6j 

Requirements (i) and (ii) are needed for algorithmic convenience only, and 
do not arise from any particular argument in § |4} Thus, in principle, this 
convenience could be traded for improvements in other areas — for example, 
for better condition numbers for the local problems, or for additional flex- 
ibility in dealing with complex geometries. However, for simplicity, in this 
paper we enforce (i) and (ii). As explained earlier (see remark [9]), we solve 
Eq. ([9]) in a least squares sense. Hence integrations over fip-' are required. It 
is thus useful to keep f2p-' as simple as possible. 

A discussion of various aspects regarding the proper definition of fip"* can 
be found in appendix [B] For instance, the requirement in item (ii) is con- 
venient only when an implicit representation of the interface is used. Fur- 
thermore, although the definition of Q^'' presented here proved robust for all 
the applications of the 4*"^ order accurate scheme (see §[6|, there are specific 
geometrical arrangements of the interface for which (ii) results in extremely 
elongated fip-'. These elongated geometries can have negative effects on the 
accuracy of the scheme. We noticed this effect in the 2^^*^ order accurate ver- 
sion of the method described in appendix [C] These issues are addressed by 
the various algorithms (of increasing complexity) presented in appendix [B] 

With the points above in mind, here (for simplicity) we define fip-' as the 
smallest rectangle that satisfies the requirements in (i), (ii), (iv), and (v) 
— then (iii) follows automatically. Hence f2p-' can be constructed using the 
following three easy steps: 

1. Find the coordinates (xminr) ^^maxr) (l/minr; 1/maxr) of the Smallest rect- 
angle satisfying condition (ii), which completely encloses the section of the 
interface F contained by the region covered by the 9-point stencil. 

2. Find the coordinates (xmin^^ x^a^^) and {ymino, yma^o) of the smallest 
rectangle satisfying condition (ii), which completely encloses all the nodes 
at which D needs to be known. 

3. Then Op-' is the smallest rectangle that encloses the two previous rectan- 
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gles. Its edges are given by 



•^max max (Xmaxp? -^maxo) 
2/min mill (?/niinp ; Z/minjj ) ; 
Z/max max (l/maxr 5 Z/max u ) 



(16a) 
(16b) 
(16c) 
(16d) 



Figure [Z] shows an example of Q^^-' defined using these specifications. 
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Figure 4: The set fJp"' for the situatfon in figure [s] 



Remark 10. Notice that for each node next to the interface we construct a 
domain ilp"'. When doing so, we allow the domains to overlap. For example, 
the domain fip-' shown in figure |4] is used to determine Cij. It should be clear 
that Q^jT^'-'^^ (used to determine Cj-ij+i), and fip'*'^'-'"^ (used to determine 



a 



each will overlap with fip"*. 



The consequence of these overlaps is that different computed values for D 
at the same node can (in fact, will) happen — depending on which domain is 
used to solve the local Cauchy problem. However, because we solve for D — 
within each fip-' — to 4*^ order accuracy, any differences that arise from this 
multiple definition of D lie within the order of accuracy of the scheme. Since 
it is convenient to keep the computations local, the values of D resulting 
from the domain f2p-', are used to evaluate the correction term Cij. X 
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Remark 11. While rare, cases where a single interface crosses the same 



stencil multiple times can occur. In § |6.3| we present such an example. A 
simple approach to deal with situations like this is as follows: First associate 
each node where the correction function is needed to a particular piece of 
interface crossing the stencil (say, the closest one). Then define one f2p"' for 
each of the individual pieces of interface crossing the stencil. 

For example, figure [sj^a) depicts a situation where the stencil is crossed 
by two pieces of the same interface (Fi and F2), with D needed at the nodes 
[i + l,j + 1), {i + l,j), + 1), and {i — l,j — 1). Then, first associate: 
(i) {i + 1, j + 1), (i + 1, j), and (i, j + 1) to Fi, and (ii) {i - 1, j - 1) to F2. 
Second, define 



1. fip-J is the smallest rectangle, parallel to the grid lines, that includes Fi 
and the nodes {i + 1, j + 1), {i + 1, j), and (i, j + 1). 

2. fip^ is the smallest rectangle, parallel to the grid lines, that includes F2 
and the node {i — 1, j — 1). 




Figure 5: Configuration where multiple fip^ are needed in the same stencil, (a) Same 
interface crossing the stencil multiple times, (b) Distinct interfaces crossing the same 
stencil. 

After the multiple Vt^ defined within a given stencil, the local Cauchy 
problem is solved for each Vfy^ separately. For example, in the case shown 
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in figure 5|a), the solution for D inside f2p-^ is done completely independent 
of the solution for D inside f2p^. The decoupling between multiple crossings 
renders the CFM flexible and robust enough to handle complex geometries 
without any special algorithmic considerations. X 



Remark 12. When multiple distinct interfaces are involved, a single stencil 
can be crossed by different interfaces - e.g.: see § |6.5| and § |6.6[ This sit- 



uation is similar to the one described in remark 11, but with an additional 



complication: there may occur distinct domain regions that are not sepa- 
rated by an interface, but rather by a third (or more) regions between them. 
An example is shown in figure |5]^b), where ri_2 and r2-3 are not part of the 
same interface. Here ri_2 is the interface between f^i and Q2, while r2-3 
is the interface between Q2 and ^3. There is no interface ri_3 separating 
Qi from Q3, hence no jump conditions between these regions are provided. 
Nonetheless, -Di_3 = (^3 — ui) is needed at (z + 1, j + 1). 

Situations such as these can be easily handled by noticing that we can 
distinguish between primary {e.g. -Di_2 and -D2-3) and secondary correction 
functions, which can be written in terms of the primary functions 
{e.g. -Di_3 = Di_2 + -D2-3) and need not be computed directly. Hence we 



can proceed exactly as in remark 11, except that we have to make sure that 
the intersections of the regions where the primary correction functions are 
computed include the nodes where the secondary correction functions are 
needed. For example, in the particular case in figure |5|(b), we define 



1. fip"'_ is the smallest rectangle, parallel to the grid lines, that includes 
ri_2 and the nodes {i + 1, j + 1), {i + 1, j), and {i,j + 1). 

^* ^ '2 — 6 

T2-2, and the node {i + 1, j + 1). 



^r2_3 smallest rectangle, parallel to the grid lines, that includes 



5.4. Solution of the Local Cauchy Problem. 

Since we use a 4*"^ order accurate discretization of the Poisson problem, we 
need to find D with 4*^ order errors (or better) to keep the overall accuracy of 
the scheme — see § 5.7 Hence we approximate D using cubic Hermite splines 

which guarantees 4**^ order accuracy 



(bicubic interpolants in 2D), which guarantees 4**^ order accuracy — see |28] 
Note also that, even though the example scheme developed here is for 2D, 
this representation can be easily extended to any number of dimensions. 
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Given a choice of basis functions]^ we solve the local Cauchy problem 
defined in Eq. ([9]) in a least squares sense, using a minimization procedure. 
Since we do not have boundary conditions, but interface conditions, we must 
resort to a minimization functional that is different from the standard one 
associated with the Poisson equation. Thus we impose the Cauchy interface 
conditions by using a penalization method. The functional to be minimized 
is then 

Jp = (£M)3 f [^^Dix)-fn{x)YdV 

+ cp [D (f) - a {x)f dS (17) 

+ cp{il:^y [ [D^ix)-b{x)fdS, 

where cp > is the penalization coefficient used to enforce the interface 
conditions, and i^:^ > is a characteristic length associated with fip-' — we 
used the shortest side length. Clearly Jp is a quadratic functional whose 
minimum (zero) occurs at the solution to Eq. 

In order to compute D in the domain fip"*, its bicubic representation is 
substituted into the formula above for Jp, with the integrals approximated 
by Gaussian quadratures — in this paper we used six quadrature points for 
the ID line integrals, and 36 points for the 2D area integrals. The resulting 
discrete problem is then minimized. Because the bicubic representation for D 
involves 12 basis polynomials, the minimization problem produces a 12 x 12 
(self-adjoint) linear system. 

Remark 13. We explored the option of enforcing the interface conditions 
using Lagrange multipliers. While this second approach yields good results, 
our experience shows that the penalization method is better. X 



Remark 14. The scaling using t^^ in Eq. (17) is so that all the three terms 
in the definition of Jp behave in the same fashion as the size of Q^'' changes 
with (i, j), or when the computational grid is refined]^ This follows because 



•^The basis functions that we use for the bicubic interpolation can be found in ap- 
pendix |Aj 

"^The scaHng also follows from dimensional consistency. 
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we expect that 





-f = o{il) 


D 


-a = 0{it) 


Dn 





Hence each of the three terms in Eq. (17) should be O 



Remark 15. Once all the terms in Eq. (17) are guaranteed to scale the same 
way with the size of f2p-', the penalization coefficient cp should be selected 
so that the three terms have (roughly) the same size for the numerical solu- 
tion (they will, of course, not vanish). In principle, cp could be determined 
from knowledge of the fourth order derivatives of the solution, which control 
the error in the numerical solution. This approach does not appear to be 
practical. A simpler method is based on the observation that cp should not 
depend on the grid size (at least to leading order, and we do not need better 
than this). Hence it can be determined empirically from a low resolution 
calculation. In the examples in this paper we found that cp ~ 50 produced 
good results. X 



Remark 16. A more general version of Jp in Eq. ( 17) would involve different 
penalization coefficients for the two line integrals, as well as the possibility of 
these coefficients having a dependence on the position along F of f2p"' . These 
modifications could be useful in cases where the solution to the Poisson prob- 
lem has large variations — e.g. a very irregular interface F, or a complicated 
forcing /. X 

5.5. Computational Cost. 

We can now infer something about the cost of the present scheme. To 
start with, let us denote the number of nodes in the x and y directions by 

Nx = ^ + l, Ny = ^ + 1, (18) 

lly 

assuming a 1 by 1 computational square. Hence, the total number of degrees 
of freedom is M = N^Ny. Furthermore, the number of nodes adjacent to the 
interface is (9(M^/^), since the interface is a ID entity. 

The discretization of Poisson's equation results in a M x M linear system. 
Furthermore, the present method produces changes only on the RHS of the 
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equations. Thus, the basic cost of inverting the system is unchanged from 
that of inverting the system resuhing from a problem without an interface F. 
Namely: it varies from 0{M) to 0{M^), depending on the solution method. 

Let us now consider the computational cost added by the modifications 
to the RHS. As presented above, for each node adjacent to the interface, 
we must construct fi^-', compute the integrals that define the local 12 x 12 
linear system, and invert it. The cost associated with these tasks is constant: 
it does not vary from node to node, and it does not change with the size 
of the mesh. Consequently the resulting additional cost is a constant times 
the number of nodes adjacent to the interface. Hence it scales as M^/^. 
Because of the (relatively large) coefficient of proportionality, for small M 
this additional cost can be comparable to the cost of inverting the Poisson 
problem. Obviously, this extra cost becomes less significant as M increases. 

5.6. Interface Representation. 

As far as the CFM is concerned, the framework needed to solve the local 
Cauchy problems is entirely described above. However, there is an impor- 
tant issue that deserves attention: the representation of the interface. This 
question is independent of the CFM. Many approaches are possible, and the 
optimal choice is geometry dependent. The discussion below is meant to shed 
some light on this issue, and motivate the solution we have adopted. 

In the present work, generally we proceed assuming that the interface 
is not known exactly - since this is what frequently happens. The only 
exceptions to this are the examples in § |6.5 and § 6^, which involve two 



distinct (circular) interfaces touching at a point. In the generic setting, in 
addition to a proper representation of the interfaces, one needs to be able 
to identify the distinct interfaces, regions in between, contact points, as well 
as distinguish between a single interface crossing the same stencil multiple 
times and multiple distinct interfaces crossing one stencil. While the CFM 
algorithm is capable of dealing with these situations once they have been 



identified {e.g. see remarks 11 and 12), the development of an algorithm with 
the capability to detect such generic geometries is beyond the scope of this 
paper, and a (hard) problem in interface representation. For these reasons. 



in the examples in § |6.5| and § |6.6| we use an explicit exact representation of 
the interface. 

To guarantee the accuracy of the solution for D, the interface conditions 



must be applied with the appropriate accuracy — see § 5.7 Since these 



conditions are imposed on the interface F, the location of F must be known 
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with the same order of accuracy desired for D. In the particular case of the 
^th qj-^Iqj- implementation of the CFM algorithm in this paper, this means 4*^ 
order accuracy. For this reason, we adopted the gradient-augmented level 
set (GA-LS) approach, as introduced in J28j. This method allows a simple 
and completely local 4**^ order accurate representation of the interface, using 
Hermite cubics defined everywhere in the domain. The approach also allows 
the computation of normal vectors in a straightforward and accurate fashion. 

We point out that the GA-LS method is not the only option for an implicit 
^th representation of the interface. For example, a regular level set 

method [27], combined with a high-order interpolation scheme, could be used 
as well. Here we adopted the GA-LS approach because of the algorithmic 
coherence that results from representing both the level set, and the correction 
functions, using the same bicubic polynomial base. 

5. 7. Error analysis. 



A naive reading of the discretized system in Eq. (14) suggests that, in 
order to obtain a fourth order accurate solution u, we need to compute the 
CFM correction terms Cij with fourth order accuracy. Thus, from Eq. (15), 



it would follow that we need to know the correction function D with sixth 
order accuracy! This is, however, not correct, as explained below. 

Since we need to compute the correction function D only at grid-points 
an 0{h) distance away from F, it should be clear that errors in the D^ j 
are equivalent to errors in a and b of the same order. But errors in a and b 
produce errors of the same order in m — see Eq. Q and Eq. ^. Hence, if we 
desire a fourth order accurate solution u, we need to compute the correction 
terms Di j with fourth order accuracy only. This argument is confirmed by 
the convergence plots in figures [7], |9| and 11 



5.8. Computation of gradients. 

Some applications require not only the solution to the Poisson problem, 
but also its gradient. Hence, in § [6| we include plots characterizing the 
behavior of the errors in the gradients of the solutions. A key question is 
then: how are these gradients computed? 

To compute the gradients near the interface, the correction function can 
be used to extend the solution across the interface, so that a standard stencil 
can be used. However, for this to work, it is important to discretize the 
gradient operator using the same nodes that are part of the 9-point stencil 
— so that the same correction functions obtained while solving the Poisson 
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equation can be used. Hence we discretize the gradient operator with a 
procedure similar to the one used to obtain the 9-point stencil. Specifically, 
we use the following 4*^^ order accurate discretization: 



where 



hi 



dxxdyUij 



dyydyUiJ 



ify)i,j 



dx'^i, 



— Ui-l,j 

2hx ' 
2^, ' 



(19) 
(20) 

(21) 
(22) 



and dxx and dyy are defined by (12) and (13), respectively. The terms {fx)i,j 
and {fy)ij may be given analytically (if known), or computed using appro- 
priate second order accurate discretizations. 

This discretization is 4*^ order accurate. However, since the error in the 
correction function is (generally) not smooth, the resulting gradient will be 
less than 4*^^ order accurate (worse case scenario is 3'''^ order accurate) next 
to the interface. 



6. Results. 

6.1. General Comments. 

In this section we present five examples of computations in 2D using the 
algorithm introduced in § [5] We solve the Poisson problem in the unit square 
[0, 1] X [0, 1] for five different configurations. Each example is defined below 
in terms of the problem parameters (source term /, and jump conditions 
across F), the representation of the interface(s) — either exact or using a 
level set function, and the exact solution (needed to evaluate the errors in 
the convergence plots). Notice that 



1. As explained in § |5.6 in examples 1 through 3 we represent the inter- 



face(s) using the GA-LS method. Hence, the interface is defined by a 
level set function 0, with gradient V(/> = {4>x,4'y) — both of which are 
carried within the GA-LS framework 1281. 
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2. Below the level set is described via an analytic formula. In examples 1 
through 3 this formula is converted into the GA-LS representation for 
the level set, before it is fed into the code. Only this representation is 
used for the actual computations. This is done so as to test the code's 
performance under generic conditions - where the interface F would be 
known via a level set representation only. 

3. Within the GA-LS framework we can, easily and accurately, compute 
the vectors normal to the interface — anywhere in the domain. Hence, 
it is convenient to write the jump in the normal derivative, [un\i^, in 
terms of the jump in the gradient of u dotted with the normal to the 
interface n = [n^, Uy). 

The last two examples involve touching circular interfaces and were devised 
to demonstrate the robustness of the CFM in the presence of interfaces that 
are very close together. In these last two examples, for the reasons discussed 
in § 5.6, we decided to use an exact representation of the circular interfaces. 



6.2. Example 1. 

• Problem parameters: 



(x, y) = — 27r^ sin(7ra;) sin(7ry), 
/~ (x, y) = — 27r^ sin(7rx) sin^Tcy), 



[u]r = sin(7rx) exp{ny), 
[un]r = [cos(7ra;) exp(7r?/) + sin(7rx) exp(7r?/) . 

• Level set defining the interface: (x, = (x — xq)^ + (y — l/o)^ — f'^i 
where xq = 0.5, yo = 0.5, and ro = 0.1. 

• Exact solution: 

(x, = sin(7rx) sin(7r?/), 
u~ (x, = sin(7rx)[sin(7r?/) — exp(7ry)]. 

Figure [6] shows the numerical solution with a fine grid (193 x 193 nodes). 
The discontinuity is captured very sharply, and it causes no oscillations in 
the solution. In addition, figure [7] shows the behavior of the error of the 
solution and its gradient in the and Loo norms. As expected, the solution 
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presents 4* order convergence as the grid is refined. Moreover, the gradient 
converges to 3^'^ order in the Loo norm and to 4*"^ order in the L2 norm, which 
is a reflection of the fact that the error in the solution is not smooth in a 
narrow region close to the interface only. 




Figure 6: Example 1 - numerical solution with 193 x 193 nodes. 




(a) Solution. (b) Gradient. 



Figure 7: Example 1 - Error behavior of the solution and its gradient in the L2 and L, 
norms. 
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6.3. Example 2. 

• Problem parameters: 



/+ (x, y) = 0, 
/" {x, y) = 0, 

[u]t = -exp(x) cos(?/), 
[un]r = ~ exp(x) cos(?/) + exp(x) sin(y) n^. 



Level set defining the interface: (x, = {x — xqY + {y — xq)^ — r'^{9), 

X — Xq 



where r{6) = ro + e sin(5^^), 6{x, y) = arctan ( - — — \ , xq = 0.5, 



= 0.5, r = 0.25, and e = 0.05. 
• Exact solution: 

(x, = 0, 
u~ (x, = exp(x) cos(y). 

Figure [s] shows the numerical solution with a fine grid (193 x 193 nodes). 
Once again, the overall quality of the solution is very satisfactory. Figure [9] 
shows the behavior of the error of the solution and its gradient in the L2 and 
Loo norms. Again, the solution converges to 4**^ order, while the gradient 
converges to 3'''^ order in the L^o norm and close to 4*^^ order in the L2 norm. 
However, unlike what happens in example 1, small wiggles are observed in the 
error plots. This behavior can be explained in terms of the construction of the 
sets ilp-' — see § Hj The approach used to construct ilp"' is highly dependent 
on the way in which the grid points are placed relative to the interface. 
Thus, as the grid is refined, the arrangement of the fl^'' can vary quite a lot 
— specially for a "complicated" interface such as the one in this example. 
What this means is that, while one can guarantee that the correction function 
D is obtained with 4*^ order precision, the proportionality coefficient is not 
constant — it may vary a little from grid to grid. This variation is responsible 
for the small oscillations observed in the convergence plot. Nevertheless, 
despite these oscillations, the overall convergence is clearly 4*^^ order. 
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Figure 9: Example 2 - Error behavior of the solution and its gradient in the L2 and L, 
norms. 



6.4- Example 3. 

• Problem parameters: 



(x, y) = exp(a;) [2 + + 2 sin(?/) + 4x sin 
/- (x, y) = 40, 
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[u]y = exp(x) [x^ sin(?/) + y^] — 10 (x^ + y"^) , 
= {exp(a;) [(x^ + 2a;) sin(?/) + y^] — 20a;} n. 



+ {exp(x) [x^ cos(|/) + 2yj — 20?/} 



Level set defining the interface: 

(x, = [(x - xi)2 + (l/ - l/i)^ - r^] [(x - X2Y + {y- - r^], wliere 
xi = yi = 0.25, ri = 0.15, X2 = 1/2 = 0.75, and r2 = 0.1. 

Exact solution: 

(x, = exp(x) [x^ sin(?/) + y^] , 
u~ (x, ?/) = 10 (x^ + y^) . 



Figure 10 shows the numerical solution with a fine grid (193 x 193 nodes). 
In this example, there are two circular interfaces in the solution domain. The 
two regions inside the circles make while the remainder of the domain 
is Q"*". This example shows that the method is general enough to deal with 



multiple interfaces, keeping the same quality in the solution. Figure [TT| shows 
that the solution converges to 4*^ order in both Loo and L2 norms, while the 
gradient converges to 3'''^ order in the Lqo norm and close to 4**^ order in the 
L2 norm. 

6.5. Example 4- 

• Problem parameters: 

/i (x, y) = -27r^ sin(7rx) sin(7ry), 

/2 (x, y) = exp(x) [2 + 1/2 + 2 sin(y) + 4x sin(y)] , 

/s (x, y) = -277^ sin(7rx) sin(7ry), 

Mri-2 = exp(x) [x^ sin(?/) + y^] — sin(7rx) sin(7ry) — 5, 

= {exp(x) [(x^ + 2x) sm{y) + y^] — 7rcos(7rx) sm(ny)}nx 
- {exp(x) [x^ cos(?/) + 2?/] — vr sin(7rx) cos(7ry)}? 



■u]r2_3 — sin(7rx)[sin(7ry) — exp(7ry)] — exp(x) [x^ sm{y) + y^] , 

= {vr cos(7rx)[sin(7ry) — exp^iry)] — exp(x) [(x^ + 2x) sm{y) + y^] }nx 
+ {7rsin(7rx)[cos(7r?/) — exp(7ry)] — exp(x) [x^ cos(?/) + 2?/] }ny. 
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Figure 10: Example 3 - numerical solution with 193 x 193 nodes. 




(a) Solution. (b) Gradient. 



Figure 11: Example 3 - Error behavior ol the solution and its gradient in the L2 and L, 
norms. 

• Interface (exact representation): 

— Region 1: inside of the big circle. 

— Region 2: outer region. 

— Region 3: inside of the small circle. 
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— Interface 1-2 (Big circle): 



tb = 0.3, 
Xob = 0.5, 
Vob = 0.5, 

Interface 2-3 (Small circle): 
rs = 0.3, 

Xqs = + rBCOs(7r/e^) - rgcos (7r(l/e^ + 1)) , 
yos =yoB +rijsin(7r/e^) - rs sin (7r(l/e^ + 1)) . 



Exact solution: 



ui (x, y) = sin(7ra;) sin^ny) + 5, 
U2 (x, y) = exp(x) [x^ sin(y) + y"^] , 
Ms (x, y) = sin(7rx)[sin(7r?/) — expi^ny)]. 



Figure 12 shows the numerical solution with a fine grid (193 x 193 nodes). 
In this example the big circle is centered within the square integration domain 
and the small circle is external to it, with a common point of tangency. The 
point of contact is placed along the boundary of the big circle at the polar 
angle 9 = vr/e^ — use the center of the big circle as the polar coordinates' 
origin. These choices guarantee that, as the grid is refined, a wide variety 
of configurations involving two distinct interfaces crossing the same stencil 
occurs in a neighborhood of the contact point. In particular, the selection 
of the angle 9 is so that no special ahgnments of the grid with the local 



geometry near the contact point can happen. Figure [13] shows the behavior 
of the error as /i — )■ in the L2 and Loo norms. Once again we observe 4*^ 
order convergence (with small superimposed oscillations) for the solution. 
Moreover, the gradient converges to 3'''^ order in the L^o norm and close to 
^th Qj-f^Qj- ^]2g norm. This example shows that the CFM is robust even 
in situations where distinct interfaces can get arbitrarily close (tangent at a 
point). 



30 




1 



Figure 12: Example 4 - numerical solution with 193 x 193 nodes. 




h h 

(a) Solution. (b) Gradient. 

Figure 13: Example 4 - Error behavior ol the solution and its gradient in the L2 and L, 
norms. 

6.6. Example 5. 

• Problem parameters: 

/i (x, y) = -27r^ sin(7rx) sin(7r?/), 
/a (x, y) = -27r^ sin(7rx) sin(7r?/), 
/s (x, y) = exp(x) [2 + + 2 sm{y) + Ax sm{y)] , 
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Mri_2 = - [sin(7rx) exp{ny) + 5] , 
^"]ri_2 ~ [cos(7ra;) exp{ny)nx + sin(7rx) exp(7r?/)R 



y} ' 



Nr2_3 = exp(x) [x^ sin(y) + y^] - sin(7ra;)[sin(7r?/) - exp(7r?/)], 
[Mn]r2_3 = {exp(x) [(x^ + 2x) sin(y) + y^] - tt cos(7rx)[sin(7r?/) - exp{7iy)]}n^ 
+ {exp(x) [x^ cos(y) + 2?/] — 7rsin(7rx)[cos(7ry) — exp(iTy)]}ny. 



Interface (exact representation): 



Region 1 
Region 2 
Region 3 



inside small circle, 
region between circles, 
outer region. 



Interface 2-3 (Big circle): 



0.3, 
0.5, 
0.5, 



Interface 1-2 (Small circle): 



rs = 0.3, 

Xqb + (i^B - rs) cos(7r/e^). 



yos 



Vob + i^B -r5)sin(7r/e^ 



Exact solution: 



Ui (x, y) = sin(7rx) sin^Tcy) + 5, 

U2 (x, y) = sin(7rx)[sin(7r?/) - exp(7r?/)], 

M3 (x, y) = exp(x) [x^ sm{y) + y"^] . 



This example complements the example in § |6.5[ the sole difference being 
that here the small circle is inside the big circle. The results are entirely 



similar. Figure 14 shows the numerical solution with a fine grid (193 x 193 



nodes), while figure 13 shows the behavior of the error in the L2 and 
norms. 
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Figure 14: Example 5 - numerical solution with 193 x 193 nodes. 




(a) Solution. (b) Gradient. 



Figure 15: Example 5 - Error behavior of the solution and its gradient in the L2 and L, 
norms. 



7. Conclusions. 

In this paper we have introduced the Correction Function Method (CFM), 
which can, in principle, be used to obtain (arbitrary) high-order accurate 
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solutions to constant coefficients Poisson problems with interface jump con- 
ditions. This method is based on extending the correction terms idea of 
the Ghost Fluid Method (GFM) to that of a correction function, defined in a 
(narrow) band enclosing the interface. This function is the solution of a PDE 
problem, which can be solved (at least in principle) to any desired order of 
accuracy. Furthermore, hke the GFM, the GFM allows the use of standard 
Poisson solvers. This feature follows from the fact that the interface jump 
conditions modify (via the correction function) only the right-hand-side of 
the discretized linear system of equations used in standard linear solvers for 
the Poisson equation. 

As an example application, the new method was used to create a 4*** or- 
der accurate scheme to solve the constant coefficients Poisson equation with 
interface jump conditions in 2D. In this scheme, the domain of definition of 
the correction function is split into many grid size rectangular patches. In 
each patch the function is represented in terms of a bicubic (with 12 free 
parameters) , and the solution is obtained by minimizing an appropriate (dis- 
cretized) quadratic functional. The correction function is thus pre-computed, 
and then is used to modify (in the standard way of the GFM) the right hand 
side of the Poisson linear system, incorporating the jump conditions into 
the Poisson solver. We used the standard 4**^ order accurate 9-point stencil 
discretization of the Laplace operator, to thus obtain a 4*^ order accurate 
method. 

Examples were computed, showing the developed scheme to be robust, 
accurate, and able to capture discontinuities sharply, without creating spuri- 
ous oscillations. Furthermore, the scheme is cost effective. First, because it 
allows the use of standard "black-box" Poisson solvers, which are normally 
tuned to be extremely efficient. Second, because the additional costs of solv- 
ing for the correction function scale linearly with the mesh spacing, which 
means that they become relatively small for large systems. 

Finally, we point out that the present method cannot be applied to all 
the interesting situations where a Poisson problem must be solved with jump 
discontinuities across an interface. Let us begin by displaying a very gen- 
eral Poisson problem with jump discontinuities at an interface. Specifically, 
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where we use the notation in § [2] and 

[31 The brackets indicate jumps across the interface, for example: 

[a m]p = (a"*" u+) (x) — (a~ u~) (x) for x eT. 

[32 The subscript n indicates the derivative in the direction of n, the unit 
normal to the interface F pointing towards 

[33 > and are smooth functions of x, defined in the union of 
and some finite width band enclosing the interface F. 

[34 > and /~ are smooth functions of x, defined in the union of Q~ 
and some finite width band enclosing the interface F. 

[35 > 0, 7^ > 0, and are smooth functions, defined on some finite 
width band enclosing the interface F. 



[36 The Dirichlet boundary conditions in (23e) could be replaced any other 
standard set of boundary conditions on dQ. 

[37 As usual, the degree of smoothness of the various data functions in- 
volved determines how high an order an algorithm can be obtained. 

Assume now that 7^ = a^, 77^ = ca^ — where c is a constant, and that 



(24) 
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applies in the band enclosing T where all the functions are defined]^ In this 
case the methods introduced in this paper can be used to deal with the 



problem in (23) with minimal alterations. The main ideas carry through, as 
follows 

[71a We assume that both m+ and u~ can be extended across F, so that they 
are defined in some band enclosing the interface. 

[7lb We define the correction function, in the band enclosing F where the 
and the exist, hj D = a^{x) {x) — a~{x) u~ [x). 

[71c We notice that the correction function satisfies the elliptic Cauchy prob- 
lem 

V^D + A- + = —/+-—/-, (25) 
with D = a and D„ = b — ca aX the interface F. 

[7ld We notice that the GFM correction terms can be written with knowl- 
edge of D. 



Unfortunately, the conditions in (24) exclude some interesting physical phe- 
nomena. In particular, in two-phase flows the case where and 7^ are 
constants (but distinct) and 77^ = arises. We are currently investigating 
ways to circumvent these limitations, so as to extend our method to problems 
involving a wider range of physical phenomena. 

A. Bicubic interpolation. 

Bicubic interpolation is similar to bilinear interpolation, and can also be 
used to represent a function in a rectangular domain. However, whereas 
bilinear interpolation requires one piece of information per vertex of the 
rectangular domain, bicubic interpolation requires 4 pieces of information: 
function value, function gradient, and first mixed derivative {i.e. fxy)- For 
completeness, the relevant formulas for bicubic interpolation are presented 
below. 

We use the classical multi-index notation, as in Ref. [2B] • Thus, we repre- 
sent the 4 vertices of the domain using the vector index v G {0, 1}^. Namely, 



Note that ( 24 ) implies that /3+ is a multiple of /3 
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the 4 vertices are = (x^ + fi Axi, X2 + V2AX2), where (x?, Xj) are the 
coordinates of the left-bottom vertex and Axj is the length of the domain 
in the Xj direction. Furthermore, given a scalar function (j), the 4 pieces of 
information needed per vertex are given by 

4 = dUix,), (A.l) 

where both v, a E {0, 1}^ and 

Then the 16 polynomials that constitute the standard basis for the bicubic 
interpolation can be written in the compact form 

2 



w^i=n<(^o, (A.3) 

1=1 

where Xj = , and is the cubic polynomial 

{/(x) for v = and a = 0, 

f (1 — x) for t> = 1 and a = 0, ,s 

( \ f n A 1 A.4 

(7(xj for 1; = U and a = 1, 

—g{l — x) for = 1 and a = 1, 

where /(x) = 1 — 3 x^ + 2 x^ and g{x) = x (1 — x)^. 

Finally, the bicubic interpolation of a scalar function is given by the 
following linear combination of the basis functions: 

nix)= 

v,de{o,i}^ 

As defined above (standard bicubic interpolation), 16 parameters are needed 
to determine the bicubic. However, in Ref. [28] a method ("cell-based ap- 
proach") is introduced, that reduces the number of degrees of freedom to 
12, without compromising accuracy. This method uses information from the 
first derivatives to obtain approximate formulae for the mixed derivatives. 
In the present work, we adopt this cell-based approach. 
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B. Issues affecting the construction of fip-'. 
B.l. Overview. 

As discussed in § [sj tlie CFM is based on local solutions to the PDE ^ 
in sub-regions of Qr ^ which we call fl^'' ■ However, there is a certain degree 
of arbitrariness in how fip"' is defined. Here we discuss several factors that 
must be consider when solving equation ([9]), and how they influence the 
choice of fip-'. We also present four distinct approaches to constructing fip-', 
of increasing level of robustness (and, unfortunately, complexity). 

The only requirements on f2p-' that the discussion in § [4] imposes are 

• f2p"' should be small, since the local problems' condition numbers in- 
crease exponentially with distance from F — see remarks [5] and |6| 

• fip-' should contain all the nodes where the correction function D is 
needed. 

In addition, practical algorithmic considerations further restrict the definition 
of fip-', as explained below. 

First, we solve for D in a weak fashion, by locally minimizing a discrete 



version of the functional Jp defined in equation (17). This procedure involves 
integrations over fip-'. Thus, it is useful if fip"' has an elementary geometrical 
shape, so that simple quadrature rules can be applied to evaluate the inte- 
grals. Second, if f2p"' is a rectangle, we can use a (high-order) bicubic (see 
appendix |a| to represent D in f2p"'. Hence we restrict fl^^'^ to be a rectangl^ 
Third, another consideration when constructing Qp-' is how the interface is 
represented. In principle the solution to the PDE in ^ depends on the infor- 
mation given along the interface only, and it is independent of the underlying 
grid. Nevertheless, consider an interface described implicitly by a level set 
function, known only in terms of its values (and perhaps derivatives) at the 
grid points. It is then convenient if the portions of the interface contained 
within f2p"' can be easily described in terms of the level set function dis- 
cretization — e.g. in terms of a pre-determined set of grid cells, such as the 



cells that define the discretization stencil. The approaches in § B.2 through 



B.4 are based on this premise. 



^Clearly, other simple geometrical shapes, with other types approximations for D, 
should be possible — though we have not investigated them. 
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Although the prior paragraph's strategy makes for an easier implemen- 
tation, it also ties Q^'' underlying grid, whereas it should depend only 
on the interface geometry. Hence it results in definitions for Q^'' that cannot 
trac k the interface optimally. For this reason, we developed the approach 
in § B.5, which allows fl^-' to freely adapt to the local interface geometry, 
regardless of the underlying grid. The idea is to first identify a piece of the 
interface based on a pre-determined set of grid cells, and then use this infor- 
mation to construct an optimal ft^ . This yields a somewhat intricate, but 
very robust definition for fip-'. 

Finally, note that explicit representations of the interface are not con- 
strained by the underlying grid. Moreover, information on the interface ge- 
ometry is readily available anywhere along the interface. Hence, in this case, 
an optimal fip-' can be constructed without the need to identify a piece of 
the interface in terms of a pre-determined set of grid cells. This fact makes 
the approach § |B.5| straightforward with explicit representations of the inter- 



face. By contrast, the less robust approaches in § B.2 through § B.4 become 
more involved in this context, because they requires the additional work of 
constraining the explicit representation to the underlying grid. 

Obviously, the algorithms presented here (§ B.2 through § B.5) represent 



only a few of the possible ways in which fip-' can be defined. Nevertheless, 
these approaches serve as practical examples of how different factors must 
be balanced to design robust schemes. 

B.2. Naive Grid-Aligned Stencil-Centered Approach. 

In this approach, we fit fip-' to the underlying grid by defining it as the 



2hx X 2hy box that covers the 9-point stencil. Figure B.l shows two examples 



This approach is very appealing because of its simplicity, but it has se- 
rious flaws and we do not recommend it. The reason is that the piece of 
the interface contained within i^p"' can become arbitrarily small - see fig- 
ure B.l[ b). Then the arguments that make the local Cauchy problem well 
posed no longer apply — see remarks [5} [6} and |8| In essence, the biggest 
frequency encoded in the interface, fcmax ~ l/length(r/f2p-'), can become ar- 
bitrarily large — while the characteristic length of fl^'' remains 0{h). As a 
consequence, the condition number for the local Cauchy problem can become 
arbitrarily large. We describe this approach here merely as an example of 
the problems that can arise from too simplistic a definition of f2p-'. 
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(a) Well-posed. (b) Ill-posed. 

Figure B.l: flp-' as defined by the naive grid-aligned stencil-centered approach. 



B.3. Compact Grid-Aligned Stencil-Centered Approach. 



This is the approach described in detail in § 5.3 In summary: f2p"' is 
defined as the smallest rectangle that 

(i) Is aligned with the grid. 

(ii) Includes the piece of the interface contained within the stencil. 

(iii) Includes all the nodes where D is needed. 



Figure B.2 shows three examples of this definition. As it should be clear 



from this figure, a key consequence of (i-iii) is that the piece of interface 
contained within fip-' is always close to its diagonal — hence it is never too 
small relative to the size of Q^'^ . Consequently, this approach is considerably 



more robust than the one in § |B.2[ In fact, we successfully employed it for 
all the examples using the 4*"^ order accurate scheme — see § [g] 

Unfortunately, the requirements (i-ii) in this approach tie f2p-' to the grid 
and the stencil. As mentioned earlier, these constraints may lead to an Q^^-' 
which is not the best fit to the geometry of the interface. Figure B.2K c) 



depicts a situation where this strategy yields relatively poor results. This 
happens when there is an almost perfect alignment of the interface with the 
grid, which can result in an excessively elongated f2p"' — in the worse case 
scenario, this set could reduce to a line. Although the local Cauchy problem 
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(c) Elongated. 



Figure B.2: rip-* as defined by the compact grid-aligned stencil-centered approach. 



remains well conditioned, the elongated sets can interfere with the process 
we use to solve the equation for Z) in ([9]). Essentially, the representation of a 
function by a bicubic (or a modified bilinear in the case of the 2°'^ accurate 
scheme in appendix [C]) becomes an ill-defined problem as the aspect ratio of 
the rectangle fip"' vanishes (however, see the next paragraph). In the authors' 
experience, when this sort of alignment happens the solution remains valid 
and the errors are still relatively small — but the convergence rate may be 
affected if the bad alignment persists over several grid refinements. 

We note that we observed the difficulties described in the paragraph above 
with the 2^^*^ accurate scheme in appendix [c] only. We attribute this to the 
fact that the bicubics result in a much better enforcement of the PDE ^ 
than the modified bilinears. The latter are mostly determined by the interface 
conditions — see remark C.l Hence, the PDE ^ provides a much stronger 
control over the bicubic parameters, making this interpolation more robust 
in elongated sets. 

Note that this issue is much less severe than the problem affecting the 
approach in 



B.2 



and could be corrected by making the Cauchy solver 
"smarter" when dealing with elongated sets. Instead we adopted the simpler 
solution of having a definition for fip-' that avoids elongated sets. The most 



robust way to do this is to abandon the requirements in (i-ii), and allow 
f2p"' to adapt to the local geometry of the interface. This is the approach 
introduced in § B.5 A simpler, compromise solution, is presented in § B.4 
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B.4- Free Stencil-Centered Approach. 

Here we present a compromise solution for avoiding elongated fip"*, which 
abandons the requirement in (i), but not in (ii) — since (ii) is convenient 
when the interface is represented implicitly. In this approach fip-' is defined 
as the smallest rectangle that 

Is aligned with the grid rotated by an angle 6^, where 6^ = 6^ — tt /4: 
and characterizes the interface alignment with respect to the grid 
{e.g. the polar angle for the tangent vector to the interface section at 
its mid-point inside the stencil). 

Includes the piece of the interface contained within the stencil, 
(iii*) Includes all the nodes where D is needed. 

Figure |B.2 shows two examples of this approach. 
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Figure B.3: ilp"' as defined by the free stencil-centered approach. 

The implementation of the present approach is very similar to that of 
The only additional work is to compute 9r 



the one in 5 B.3 



and to write 

the interface and points where D is needed in the rotated frame of reference. 
In both these approaches the diagonal of f2p"' is very close to the piece of 
interface contained within the stencil, which guarantees a well conditioned 
local problem. However, here the addition of a rotation keeps f2p"' nearly 
square, and avoids elon gated geometries. The price paid for this is t hat t he 
sets fip-' created using 



B.4 



B.3 



can be a little larger than the ones from 
with both sets including the exact same piece of interface. In such situations, 
the present approach results in a somewhat larger condition number. 
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B.5. Node- Centered Approach. 

Here we define in a fashion tliat is completely independent from 
the underlying grid and stencils. In fact, instead of associating each Q^^-' to a 
particular stencil, we define a different f2p"' for each node where the correction 
is needed — hence the name node-centered, rather than stencil-centered. As 
a consequence, whereas the prior strategies lead to multiple values of D at 



1 



the same node (one value per stencil, see remark 10), here there is a unique 
value of D at each node. 

In this approach, fip"* is defined by the following steps: 

Identify the interface in the 4 grid cells that surround a given node. This 
step can be skipped if the interface is represented explicitly. 

2. Find the point, Pq, along the interface that is closest to the node. This 
point becomes the center of There is no need to obtain Pq very 
accurately. Small errors in Pq result in small shifts in ilp"' only. 

3. Compute to, the vector tangent to the interface at Pq. This vector de- 
fines one of the diagonals of fl^'' ■ The normal vector no defines the other 
diagonal. Again, high accuracy is not needed. 

4. Then fip-' is the square with side length l-^hl + K^, centered at Pq, and 
diagonals parallel to to and fiQ — fip"* need not be aligned with the grid. 



Figure B.4 shows two examples of this approach. 
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Figure B.4: Vfy'' as defined by the node-centered approach. 
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Note that the piece of interface contained within fip-', as defined by steps 
1-4 above, is not necessarily the same found in step 1. Hence, after defining 
f2p-', we still need to identify the piece of interface that lies within it. For an 
explicit representation of the interface, this additional step is not particularly 
costly, but the same is not true for an implicit representation. 

This approach is very robust because it always creates a square fip-', with 
the interface within it close to one of the diagonals — as guaranteed by 
steps 2 and 3. Hence, the local Cauchy problem is always well conditioned. 
Furthermore, making squ are ( to avoid elongation) does not result in 

because the larger f2p-' contain an equally 



B.4 



larger condition numbers as in 
larger piece of the interface. 

Finally, we point out that the (small) oscillations observed in the con- 
vergence plots shown in § |6] occur because in these calculations we use the 



approach in § B.3 — which produces sets f2p-' that are not uniform in size, 
nor shape, along the interface. Tests we have done show that these oscilla- 
tions do not occur with the node-centered approach here, for which all the 
f2p-' are squares of the same size. Unfortunately, as pointed out earlier, the 
node-centered approach is not well suited for calculations using an interface 
represented implicitly. 

C. 2"*^ Order Accurate Scheme in 2D. 

C.l. Overview. 

In § [s] we present a 4*^^ order accurate scheme to solve the 2D Poisson 
problem with interface jump conditions, based on the correction function 
defined in § |4| However, there are many situations in which a 2°^^ order 
version of the method could be of practical relevance. Hence, in this section 
we use the general framework provided by the correction function method 
(CFM) to develop a specific example of a 2°'^ order accurate scheme in 2D. 
The basic approach is analogous to the 4*^ version presented in § [s] A few 
key points are 

(a) We discretize Poisson's equation using the standard 5-point stencil. This 
stencil is compact, which is an important requirement for a well condi- 
tioned problem for the correction function D — see remarks |4]-[6] 



(b) We approximate D using modified bilinear interpolants (defined in § C.2 ) 



each valid in a small neighborhood Vt^ ^^e interface. This guarantees 



local 2 order accuracy with only 5 interpolation parameters. Each 
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fip-' corresponds to a grid point at which the standard discretization of 
Poisson's equation involves a stencil that straddles the interface F. 

(c) The domains Q^'^ are rectangular regions, each enclosing a portion of F, 
and all the nodes where D is needed to complete the discretization of 
Poisson's equation at the (i, j)-th stencil. Each is a sub-domain of Qr- 

(d) Starting from (b) and (c), we design a local solver that provides an 
approximation to D inside each domain fip-'. 

(e) The interface F is represented using the standard level set approach — 
see [27]. This guarantees a local 2°*^ order representation of the interface, 
as required to keep the overall accuracy of the scheme. 

(f) In each fip-', we solve the PDE in ([9]) in a least squares sense — see re- 
mark |9} Namely, we seek the minimum of the positive quadratic integral 



quantity Jp in (17), which vanishes at the solution: We substitute the 
modified bilinear approximation for D into Jp, discretize the integrals 
using Gaussian quadratures, and minimize the resulting discrete Jp. 

Remark C.l. Using standard bilinear interpolants to approximate D in 
each also yields 2^^*^ order accuracy. However, the Laplacian of a standard 
bilinear interpolant vanishes. Thus, this basis cannot take full advantage of 
the fact that D is the solution to the Cauchy problem in ([9]). The modified 
bilinears, on the other hand, incorporate the average of the Laplacian into the 



formulation — see 5 C.2 In 5 C.6 we include results with both the standard 



and the modified bilinears, to demonstrate the advantages of the latter. X 



Below, in § C.2 C.5 we describe the 2°^^ order accurate scheme, and in 



C.6 we present applications of the scheme to three test cases. 



C.2. Modified Bilinear. 

The modified bilinear interpolant used here builds on the standard bi- 
linear polynomials. To define them, we use the multi-index notation — 
see appendix |A] and Ref. [2H]- Thus, we label the 4 vertices of a rect- 
angular cell using the vector index v G {0, 1}^. Then the vertices are 
Xv = (x? + f 1 Axi, X2 -\- V2 Aa;2), where ( the coordinates of the 

left-bottom vertex and Axi is the length of the domain in the Xi direction. 
The standard bilinear interpolation basis is then given by the 4 polynomials 

2 

W^^ = l[w^^{x,), (C.l) 
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where Xi 



r.0 



Axi ' 



and w"" is the hnear polynomial 



w^ix) 



1 ~ X for V = 0, 
X for V = 1. 



(C.2) 



The standard bilinear interpolation of a scalar function (j) is given (in each 
cell) by 

'Hs{x)= Yl WU{x^)- (C.3) 

v€ {0,1}2 

In the modified version, we add a quadratic term proportional to + y^, so 
that the Laplacian of the modified bilinear is no longer identically zero. The 
coefficient of the quadratic term can be written in terms of the average value 
of the Laplacian over the domain, 0. This yields the following formula 
for the modified bilinear interpolant 



Hm (x) = Hs [x) - - [w\x)w\x){Axf + w\y)w\y){/\yf] 0. (C.4) 

C.3. Standard Stencil. 

We use the standard 2^^^^ order accurate 5-point discretization of the Pois- 
son equation 



(C.5) 



where L is defined in (11). In the absence of discontinuities, (C.5) provides 



a compact 2°'^ order accurate representation of the Poisson equation. In the 
vicinity of the discontinuities at the interface P, we define an appropriate 



domain f2p-', and use it to compute the correction terms needed by (C.5) — 
as described below. 

To understand how the correction terms affect the discretization, consider 



the situation in figure C.l In this case, the node («, j) lies in fi"*" while the 



nodes (z+1, j) and (i, j + 1) are in Vt~ . Hence, to be able to use equation (C.5 ) 
we need to compute -Dj+ij and -Djj+i- 

After having solved for D where necessary (see 



modify equation (C.5) and write 



C.4 and § C.5) 



we 



(C.6) 



which differs from (C.5) by the terms Cij on the RHS only. Here the Cij 
are the CFM correction terms needed to complete the stencil across the 
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Figure C.l: The 5-point stencil next to the interface F. The dashed box shows a compact 
quadrangular region that contains the stencil. 



discontinuity at F. In the particular case illustrated by figure C.l , we have 



(C.7) 



Similar formulas apply for all the other possible arrangements of the stencil 
for the Poisson's equation, relative to the interface F. 

C.4- Definition of Vt^ ■ 



As discussed in appendix |B| the construction of Vfy-' presented in § 
may lead to elongated shapes for r2 p-', w hich can cause accuracy loses, 
mentioned earlier (see the end of 



5.3 



As 



B.3) this is not a problem for the 4*^ 
order scheme, but it can be one for the 2°*^ order one. To resolve this issue 



we could have implemented the robust construction of Vt'-' 



P given m 



B.5 



However the simpler c omp romise version in § B.4 proved sufficient. 



The approach in 



B.4 



tained "within the stenci' 



requires Vt-f to include the piece of interface con- 
." For the 9-point stencil, this naturally means 
"within the 2hx x 2hy box aligned with the grid that includes the nine points 
of the stencil." We could use the same meaning for the 5-point stencil, but 
this would not take full advantage of the 5-point stencil compactness. A 
better choice is to use the quadrilateral defined by the stencil's four extreme 
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points — i.e. the dashed box in figure C.l This choice is used only to de- 
termine the piece of interface to be included within Op-'. It does not affect 
the definition of the interface alignment angle di used by the §|B .41 approach. 



Hence the following four steps define Q^-' 

1. Find the angle 6^ between the vector tangent to the interface at the mid- 
point within the stencil, and the x-axis. Introduce the coordinate system 
p-q, resulting from rotating x-y by 6*^ = 6'r — vr/4 — see figure C.2[ a). 

2. Find the coordinates (pminr; Pmaxr) and (gminr) Q'maxr) characterizing the 
smallest p-q coordinate rectangle enclosing the section of the interface 
contained within the stencil — see figure C.2[ b). 

3. Find the coordinates (pmino, Pmaxo) and {qmmo, ?maxo) characterizing the 
smallest p-q coordinate rectangle enclosing all the nodes at which D is 
needed — see figure C.2[ b). 



fip-' is the smallest p-q coordinate rectangle enclosing the two previous 



rectangles. Its edges are characterized by 



nin min (Pminr ? Pmini5 ) ; 
Piaax 



max (Pmaxp ; Pmax£3 ) ; 



Q'min min (o'minp) Q'min/j 
9max max (O'maxr! Q'max/3 y 



(C 
(C.8b) 
(C.Sc) 
(C.8d) 



Figure C.2 shows an example of fip"' defined in this way. 



C.5. Solution of the Local Cauchy Problem. 

The remainder of the 2^^*^ order accurate scheme follows the exact same 
procedure used for the 4**^ order accurate version described in detail in § [s} In 
particular, as explained in item (f) of § C.l , we solve the local Cauchy problem 



defined by equation (|9]) in a least squares sense (using the same minimization 
procedure described in § 5.4). The only differences are that: (i) In the 2°^ 



order accurate version of the method we use 4 Gaussian quadrature points 
for the ID line integrals, and 16 points for the 2D area integrals, (ii) Because 
the modified bilinear representation for D involves 5 basis polynomials, the 
minimization problem produces a 5 x 5 (self-adjoint) linear system — instead 
of the 12 X 12 system that occurs for the 4*^order algorithm. 
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(a) Step 1. (b) Steps 2-4. 



Figure C.2: The set fJp"' for the situatfon in figure C.l 



C. 6. Results 

Here we present three examples of solutions to the 2D Poisson's equa- 
tion using the algorithm described above. Each example is defined below in 
terms of the problem parameters (source term / and jump conditions across 
the interface F), the representation of the interface, and the exact solution 
(needed to evaluate the errors in the convergence plots). Note that 

1. In all cases we represent the interface using a level set function cf). 

2. Below the level set function is described via an analytic formula. This 
formula is converted into a discrete level set representation used by the 
code. Only this discrete representation is used for the actual computa- 
tions. 

3. The level set formulation allows us to compute the vectors normal to the 
interface to 2^'^ order accuracy with a combination of finite differences 
and standard bilinear interpolation. Hence, it is convenient to write 
the jump in the normal derivative, [m^Jp, in terms of the jump in the 
gradient of u dotted with the normal to the interface n = (n^, ny). 
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We us^ example 1^ to test the 2'^'^ order scheme in a generic problem with 
a non-trivial interface geometry. In addition, in this example the correction 
function D has a non-zero Laplacian. Thus, we used this problem to compare 
the performance of the Standard Bilinear (SB) and Modified Bilinear (MB) 
interpolations as basis for the correction function. Finally, examples 2^ and 
3s here correspond to examples 1 and 3 in [3], respectively. Hence they can 
be used to compare the performance of the present 2'^'^ order scheme with 
the Immersed Interface Method (IIM), in two distinct situations. 

C.6.1. Example Is- 

• Domain: {x,y) G [0,1] x [0,1]. 

• Problem parameters: 

4, 
0, 

+ — exp(a;) cos{y), 
[2x — exp(x) cos(?/)] Tlx + [2y + exp(a;) sm{y)] Uy. 

• Level set defining the interface: (f){x, y) = ^/ {x — xq) + {y — Xo)—r(6), 

where r(6) = tq + e sin(5 6'), 6(x, y) = arctan ( - — — ) , a;o = 0.5, 

\x - XoJ 

yQ = 0.5, r = 0.25, and e = 0.05. 

• Exact solution: 

u+ (x, y) = x^ + y^, 
u~ (x, y) = exp(a:;) cos(?/). 



/ {x, y) = 
f~ {x, y) = 




Figure C.3 shows the numerical solution with a fine grid (193 x 193 nodes). 
The non-trivial contour of the interface is accurately represented and the 
discontinuity is captured very sharply. For comparison we solved this problem 
using both the SB and MB interpolations to represent the correction function. 
Although figure C.3 shows the solution obtained with the modified bilinear 
version, both versions produce small errors and are visually indistinguishable. 



subscript s is added to the example numbers of the 2'^'^order method, to avoid 
confusion with the 4"^order method examples. 
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Figure Q4](a) shows the convergence of the solution error in the L2 and 
Loo norms, for both the SB and MB versions. As expected, the overall 
behavior indicates 2°*^ order convergence, despite the small oscillations that 



are characteristic of this implementation of the method, as explained in § |6.3 
Note that the MB version produces significantly smaller errors (a factor of 
about 20) than the SB version. It also exhibits a more robustly 2"*^ order 
convergence rate. Moreover, figure C.4[ b) shows the convergence of the error 
of the gradient of the MB solution in the L2 and L^o norms. Here, the gradient 
was computed using standard 2°^^ order accurate centered differences, as in 



(21) and (22). As we can observe, the gradient converges to order in the 



Loo norm and, apparently, as h^^"^ in the L2 norm. 




Figure C.3: Example Is- 
modified bilinear version. 



Numerical solution with 193 x 193 nodes, 2"'^order scheme 
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Figure C.4: Example 1^ - Error behavior of the solution and its gradient in the L2 and 
Lao norms. 



C.6.2. Example 2^. 

• Domain: {x,y) G [—1,1] x [—1,1]. 



Problem parameters: 



t (a;, y) = 0, 
/" (a;, y) = 0, 

[M]r = log (2^/x^ + y^y 

_ xn^ + yuy 



\u. 



nJr 9,9" 
+ y^ 

Level set defining the interface: (x, y) = ^/ {x — Xq)^ + {y ~ y^Y — ro, 
where xq = 0, ?/o = 0, and ro = 0.5. 

Exact solution: 

M+ (x, ?/) = 1 + log ^2v^x2 + 1/2^ , 
u~ (x, y) = 1. 
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As mentioned earlier, this example corresponds to example 1 in [H], where 
the same problem is solved using the Immersed Interface Method. Hence, this 
provides a good opportunity to compare the CFM with the well established 
IIM. Figure C.5 shows the numerical solution with a fine grid (161 x 161 
nodes). Once again, the overall quality of the solution is very satisfactory. 
Figure C.6 shows the behavior of the error in the L2 and L^o norms and the 
solution converges to 2^^^^ order. However, in this example the gradient also 
appears to converge to 2^^^^ order. 

Figure C.6K a) also includes the convergence of the solution error in the 
Loo norm obtained with the IIM — we plot the errors listed in table 1 of [3] . 
Both methods produce similar convergence rates. However, in this example 
at least, the CFM produces slightly smaller errors — by a factor of about 
1.7. 




Figure C.5: Example 2^. Numerical solution with 161 x 161 nodes. 



53 




(a) Solution. (b) Gradient. 



Figure C.6: Example 2^ - Error behavior of the solution and its gradient in the L2 and 
Lao norms. IIM refers to results obtained with the Immersed Interface Method in [3] 
(Copyright ©1994 Society for Industrial and Applied Mathematics. Data compiled with 
permission. All rights reserved). 



C.6. 3. Example 3s. 

• Domain: {x,y) G [—1,1] x [—1,1]. 

• Problem parameters: 

{x, y) = 0, 
/" {x, y) = 0, 

[u]r = -exp{x) cos{y), 
[un]r = exp(a;) [- cos{y) + sm{y) Uy] . 



Level set defining the interface: (x, y) = \/ {x — Xq)^ + {y — VoY — fo, 
where Xq = 0, yo = 0, and tq = 0.5. 



Exact solution: 



(x, y) = 0, 
u~ (x, y) = exp(x) cos{y). 



This example corresponds to example 3 in [3] for the IIM method. Figure C.7 



shows the numerical solution with a fine grid (161 x 161 nodes) while fig- 
ure C.8 presents convergence results for the errors in the L2 and L^o norms. 
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In addition, figure C.8[ a) also includes the norm of the error obtained 
with the JIM — we plot the errors listed in table 3 of [3J In this case the IIM 
produces slightly smaller errors than the CFM — by a factor of about 2.8. 
Therefore, based on the results from examples 2^ and 3^, we may conclude 
that the IIM and the CFM produce results of comparable accuracy — each 
method generating slightly smaller errors in different cases. 




y -1 -1 



Figure C.7: Example 3s- Numerical solution with 161 x 161 nodes. 
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10 
h 



fa) Solution. 



(b) Gradient. 



Figure C.8: Example 3^ - Error behavior of the solution and its gradient in the L2 and 
Lao norms. IIM refers to results obtained with the Immersed Interface Method in [3] 
(Copyright ©1994 Society for Industrial and Applied Mathematics. Data compiled with 
permission. All rights reserved). 
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